import numpy as np


#fun= ConDfai(f,x,x0)%
#grad = dConDfai(f,x,x0)	%求f(x)的梯度
#hesse = hessen(f,x,x0)    %计算 hessen 矩阵
#P = -grad; % 选取初始下降方向
#t = (grad*(grad'))/(grad*hesse*(grad')); %预处理--第一次迭代
#x0 = x0 + t*P;
#y = x0;
#k = 1; %设定初值k=1，限制计算共轭方向次数
#while norm(grad)>=epsilon%精度判断条件
#    alpha = (grad*Q*P')/(P*Q*P'); % 计算共轭方向 al
#    P = -grad + alpha*P;%生成新的搜索方向
#    t = (grad*(grad'))/(grad*hesse*(grad'));%求出此时的步长
#    y = y + t*P;%求出新的x（k）
#    grad = dConDfai(f,x,y);     
#%     if k == n-1% k 值限制计算共轭方向次数
#%         break
#%     end
#    k = k + 1;
#end
#x0 = y;%将求出的最后极小点重新赋予x0

#syms x1 x2 x3    %x4 x5 x6 x7 x8 x9; % 初始化变量 x
#epsilon = 0.00001;   %定义精度
#x0 = [1 0 0] ; % 定义初始点
#Q = [2 1 -1; 1 4 0; -1 0 2]; % 要求 Q是对称正定矩阵                 
#b = [-1 1 2];c = 0; n = length(x0);% 求出最大维数
#x = [x1 x2 x3]; %此时维数为3
#f =0.5* x*Q*x.'+ b*x.'+ c; % 定义目标函数并计算初始值
#


def conjugatedirection(fun,gfun,hfun,Q,P,x0,eps,k,maxk):
    y = fun(x0)
    grad = gfun(x0)
    hess = hfun(x0)
    if np.lingla.norm(grad) < eps or k >= maxk:
        return x0,fun(x0),k
    alpha = (grad*Q*P.T)/(P*Q*P.T)
    P = -grad + alpha*P
    t = (grad*(grad.T))/(grad*hess*(grad.T))
    return conjugatedirection(fun,gfun,hfun,P,x0+P*t,eps,k+1,maxk)

gfunc = lambda x: np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)])

result = gradient_descent(func = lambda x:100*(x[1]-x[0])**2 + (1-x[0])**2,
                        gfunc = gfunc,
                        Q = [[2,1,-1],
                            [1,4,0],
                            [-1,0,2]] 
                        ,P = -gfunc([1,0,0]),
                        x0=[1,0,0],eps=1e-5,k=0,maxk=5000)
print(result)
